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ABSTRACT 

Popular models for describing the luminosity-density profiles of dynamically hot stellar 
systems (e.g., Jaffe, Hernquist, Dehnen) were constructed to match the deprojected 
form of de Vaucouleurs' i? 1 / 4 light-profile. However, we now know that elliptical galax- 
ies and bulges display a mass-dependent range of structural profiles. To compensate 
this, the model in Terzic & Graham was designed to closely match the deprojected 
form of Sersic i? 1 /™ light-profiles, including deprojected exponential light-profiles and 
galaxies with partially depleted cores. It is thus applicable for describing bulges in spi- 
ral galaxies, dwarf elliptical galaxies, both "power-law" and "core" elliptical galaxies, 
also dark matter halos formed from ACDM cosmological simulations. In this paper, we 
present a new family of triaxial density-potential-force triplets, which generalizes the 
spherical model reported in Terzic & Graham to three dimensions. If the (optional) 
power-law core is present, it is a 5-paramcter family, while in the absence of the core 
it reduces to 3 parameters. The isodensity contours in the new family are stratified 
on confocal ellipsoids and the potential and forces are expressed in terms of integrals 
which are easy to evaluate numerically. We provide the community with a suite of 
numerical routines for orbit integration, which feature: optimized computations of po- 
tential and forces for this family; the ability to run simulations on parallel platforms; 
and modular and easily editable design. 

Key words: galaxies: elliptical and lenticular, cD - galaxies: structure - galaxies: 
nuclei - Galaxy: bulge - galaxies: halos - dark matter 



1 INTRODUCTION 



Elliptical galaxies, bulges of spiral galaxies and dark mat- 
ter halos exhibit a range of density-profile shapes which are 
accurately described by the Sersic (1963, 1968) R 1/n model 
(e.g. Caon, Capaccioli & D'Onofrio 1993; Young & Currie 
1994; Graham et al. 1996; Graham 2001; Balcells et al. 2003; 
Merritt et al. 2006). This model is a generalization of de Vau- 
couleurs' (1948, 1959) R 1 ' 4 model which provides a decent 
fit only to galaxies with Mb ~ —21 mag (e.g. Kormendy & 
Djorgovski 1989; Graham & Guzman 2003). The popularity 
of de Vaucouleurs' model permeated into computer modeling 
of elliptical galaxies: the most popular models used in galaxy 
dynamics simulations - Jaffe (1983), Hernquist (1990) and 
Dehnen (1993, see also Tremaine et al. 1994) - were designed 
to approximate the R}^ A light profile when projected. All 
three of these density models are double power-law models 
and have the same fixed outer profile slope, declining with 
radius as r~ 4 . This means that, while undeniably useful, 
these models are limited in their ability to model realistic 



galactic structures and their evolution (Terzic & Graham 
2005). 

More sophisticated power-law models have since been 
developed to allow for more flexibility in matching depro- 
jected galaxy light profiles. A 3-parameter Dehnen double 
power- law model was generalized by Zhao (1996) into a 5- 
parameter model such that both the inner and outer power- 
law slope can be adjusted, along with the radius, density 
and sharpness of the transition region. This model was first 
introduced by Hernquist (1990, his equation 43) and has the 
same structural form as the "Nuker" model (e.g., Lauer et 
al. 1995). However, Graham et al. (2003) show that such 
double power-law models do not provide an adequate de- 
scription of profiles with logarithmic curvature, i.e. profiles 
without inner and outer power-laws which is the case for 
the luminosity-density profiles of most elliptical galaxies and 
bulges. The 3-parameter model of Prugniel & Simien (1997) 
was, however, designed to describe profiles with curvature 
and does therefore not suffer from such a problem. Their 
density model provides a good analytic approximation to 
the deprojected Sersic R 1 ^" model. 
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These expressions have proven useful for describing the 
density profiles of dark matter halos simulated in cosmo- 
logical models of hierarchical structure formation (Merritt 
et al. 2006; Graham et al. 2006). Merritt used techniques 
from nonparametric function estimation theory to extract 
the density profiles, and their derivatives, from a set of N- 
body halos, to show that the 3-parameter Prugniel-Simien 
density profile provides a better description to the data than 
other 3-parameter models currently in use, such as the gen- 
eralized NFW (Navarro, Frenk, & White 1996) model with 
arbitrary inner slope. 

Recently, Terzic & Graham (2005) developed a 5- 
parameter model which is able to unite properly the do- 
main of the black hole with the outer parts of the galaxy. It 
modifies the 3-parameter Prugniel-Simien density model to 
allow for (optional) partially depleted power-law cores, pro- 
viding excellent match to the deprojected light profiles of 
elliptical galaxies, including the ones with depleted cores for 
which the Prugniel-Simien model is inadequate. Expressions 
for the potential and forces corresponding to the spherical 
Prugniel-Simien model are reported in Terzic & Graham 
(2005). 

Double power-law models owe a good deal of their pop- 
ularity to the fact that the expressions for the potential and 
forces are relatively simple, in some cases even expressible 
in terms of analytic functions. This leads to easy implemen- 
tation of computer routines for orbital integration. In the 
trade-off between faithfulness to the physical problem (as 
measured by how accurately these density profiles match 
the observed deprojected light-profiles of elliptical galaxies) 
on one side, and numerical simplicity and computational ef- 
ficiency on the other, these models favor the latter. However, 
after more than two decades of rapid progress in computer 
technology since the first of these models was devised, this 
choice becomes more difficult to justify. Here we propose to 
bias the trade-off in the opposite direction - we devise a more 
realistic physical model at the cost of making it mathemat- 
ically more intricate and computationally more demanding. 
Our rationale is that if we are to construct realistic mod- 
els of galaxies, it is imperative that we use a model capable 
of reflecting, as accurately as possible, what is observed in 
Nature. We then estimate computational efficiency of this 
realistic model and compare it to that of traditional power- 
law models. Finally, we design and make freely available an 
optimized suite of computer programs used for orbit inte- 
gration, in hope to eliminate numerical implementation as 
an adverse factor. 

In this paper we generalize spherical profiles of Terzic 
& Graham (2005) to three dimensions by replacing radial 
with ellipsoidal symmetry in the mass density. The outline 
of the new triaxial model in the context of our earlier work 
is shown in Fig. [1] In Section [2] the new triaxial model is 
presented: we start with ellipsoidally stratified mass den- 
sity distribution and derive the corresponding potential and 
forces, including special cases when the 5-parameter model 
reduces to a 3-parameter deprojected Sersic profile with no 
power- law core in Section ^, 3l and 2-parameter power-law in 
Section 12.41 The potential and forces are derived in terms 
of integrals , which are easy to evaluate numerically. In Sec- 
tion |3] we outline the main features of the numerical suite 
used for orbital integration in the new potential, including 
code optimization and parallelization, user interface and in- 



formation analysis. In Section [4] we use the new model to 
describe the bulge at the center of the Milky Way as well as 
a pair of galaxy profiles. We also compare its speed and the 
goodness of fit to the popular Dehnen model (|Dehnen 19931 
ITremaine et al. 1994|) . Finally, we summarize and discuss 
the importance of the model presented here in Section [S] 



2 TRIAXIAL MODEL 

In this Section we introduce a family of expressions asso- 
ciated with the spatial density profiles of triaxial stellar 
systems having Sersic-like profiles with optional power-law 
cores, generalizing our earlier work on equivalent spherically 
symmetric profiles (Terzic & Graham 2005). We derive exact 
expressions for the potential and forces associated with this 
family of spatial density profiles in terms of quadratures. 
This profile subsumes two other profiles as its special cases: 
the 3-parameter Prugniel-Simien profile in the limit of the 
break radius r;, — > 0, and the 2-parameter power-law core in 
the limit of r% — > oo. 



2.1 Density 

The triaxial mass density model is a triaxial generalization 
to the spherical model of Terzic & Graham (2005) model. 
Replacing radial dependence with a dependence on an ellip- 
soidal coordinate m, defined as 



(1) 



yields a profile with equidensity surfaces stratified on confo- 
cal ellipsoids with axis ratios a : b : c 



p(m) 



if m < rt 
if m > r& 



with 



(->) 



(3) 



Parameters above are defined in the context of the spher- 
ical model and its fits to the observed light profiles 
l|Terzic fc G raham 2005 ) . The break radius rj, marks the po- 
sition of instantaneous transition from one regime to an- 
other, is the density at the break radius. 7 is the in- 
ner power-law slope, not to be confused with incomplete 
gamma function j[a, x] defined later in the text. R e is the 
(projected) effective half-light radius when r;, = 0. The pa- 
rameter n is the Sersic index and describes the curvature of 
the profile. The term b n is not a parameter, but instead a 
function of n and chosen to ensure R D contains half the (pro- 
jected) galaxy light. It is obtained by solving the equation 
r[2n] =2 j[2n,b n ], where 



T[a] = 



j[a,x\ 



T[a,x] = F[a]--y[a,x]= / e'h^dt, a>0 



(4) 
(5) 
(6) 
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Figure 1. Flow-chart outlining the work presented here. 



are the complete, incomplete and complement to the in- 
complete gamma functions, respectively. A good approxi- 
mation of b n for n > 0.5 is 2n - 1/3 + 0.009876/n (Prug- 
niel & Simien 1997; see MacArthur, Courteau, & Holtzman 
2003 for smaller values of n). Similarly to b n , p is also not 
a parameter, but a function of n, well approximated by 
p = 1.0 - 0.6097/n + 0.05563/n 2 , for 0.6 < n < 10. 



2.2 Potential and Forces 

The potential arising from an ellipsoidally stratified density 
distribution p — pirn 2 ) can be written as (Chandrasekhar 
1969, page 52, Theorem 12): 



&(x,y,z) — —nabcG 



ip(oo) — tp(m 2 ) 



v /( T + a 2)( r + b 2)( r + c 2) 



with 



m 2 (r) 



and 



ip(m 2 ) 



x 2 y 2 z 2 
+ -nr- + ■ 



a 2 + r b 2 + t c 2 + r ' 



p(rh 2 ) dm 2 . 



dr, (7) 



(8) 



(9) 



After some algebra, we find for 7 7^ 2: 
®(x,y,z) = —^—r',lH(x,y,z;a 2 ,b 2 ,c 2 ,2-^) 



+ aKR(x,y,z;a 2 ,b 2 ,c 2 ), 



when m < r b , and 



(10) 



<&(x,y,z) = ■^^—r^H I (x,y,z;a 2 ,b 2 ,c 2 ,2-j,T b ) 
2-7 



+ aK R 1 (x,y, z; a 2 ,b 2 ,c 2 ,Th.) 



ipP 2 

when m > rt, where 
a = 2-KabcGpt 

K = --^r 2 b -npR 2 b^- 2) r 
2-7 

H'(x,y,z;A,B,C,q,s) 



anpRtb" KP z> T(x,y,z;a 2 ,b 2 ,c 2 ,T b ), (11) 



n(2 



■->■*• (2) 



l/r, 



m q dr 



(12) 
,(13) 



(14) 



di{r;A,B,Cy 

S 

H(x, y, z; A, B, C, q) = H 1 \x, y, z; A, B, C, q, 0), (15) 
dr 



Rl{A > B > c > s ^Jl^T-A,B,cy 

s 

R{A,B,C) = R'{A,B,C,Q), 



(16) 
(17) 
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T(x,y,z;A,B,C,s) = 



and 



n(2-p),fc,(£) 



1 / r, 



(It 



di(r;A,B,C) 



(18) 



(19) 



di(r;A,B,C) = y/(r + A)(t + B)(t + C). 

Superscript / denotes an improper integral with a finite non- 
zero lower limit and infinity as an upper limit. R(A, B,C) 
is known as Carlson integral (Carlson 1988). r b is found by 
solving a cubic equation rt = m(T b ). 

For 7 = 2, the expressions for potential are: 



$(x,y,z) 



ar b Hi(x,y, z;a ,b , c 



(20) 



+ aK2R(a ,b , c ), 
when m < r&, and 

$>(x,y,z) = arlH[(x,y,z;a 2 ,b 2 ,c 2 ,T b ) 
+ aKiR 1 (a 2 ,b 2 , c 2 ,r b ) 

- anpR 2 b" (p ~ 2) T(x,y,z;a 2 ,b 2 ,c 2 ,T b ), (21) 
when m > r b , where 



K 2 = -rt log r b -npR 2 b^ v 2) Y 
and 

oo 

logm dr 



n{2 - p), &„ ^ 



l/n 



(22) 



H{(x,y,z\A,B,C,s) 



di(r;A,B,Cy 



(23) 
(24) 

(25) 



tf,(z, 2/, A B, C) = H((x, y, z; A, B, C, 0). 

The gravitational forces for 7 7^ 2 are given by 
2) = -ar b 1 q i H 2 (x,y,z;k 2 ,a 2 ,b 2 ,c 2 ,-j), 
when m < rb and 

F q ,{x,y,z)= - ar2q t Hi(x,y,z;k 2 ,a 2 ,b 2 ,c 2 ,~-f,r b ) 

- apq l S{x,y,z\k 2 ,a 2 ,b 2 ,c,T b ), (26) 

when m > rv fci = a, k 2 = b, ks — c and q\ = x, q 2 — y, 
q-i — z. The quadratures in the above expressions are defined 

a.s 



00 

Hi(x,y,z;D,A,B,C,q,s) = J 



m q dr 



d 2 (r;D,A,B,C)' 



(27) 



H 2 (x, y, z; D, A, B, C, q) = H' 2 {x, y, z; D, A, B, C, q, 0), (28) 



S(x,y,z;D,A,B,C,s) 



(f) 



-p —b n 

e 



dr 



d 2 (r;D,A,B,C) 



; (29) 



where 

d 2 (r;D,A,B,C) = (r + D) di(r;A,B,C). 



(30) 



For 7 = 2, no special care in derivation of forces is 
needed, so they coincide with equations (|25|l and (12611 , 

The expressions Q2), CESJ), ([23J, (J2TJ) and <[T8j) , as well 
as (I29|) when s — » 00, are improper integrals and are not very 
suitable for numerical computation. This problem is resolved 
by introducing a new variable of integration £ = (1 + t) _1//2 . 
The integrals then become proper and are given by 



H I (x,y,z;A,B,C,q,s) = 2 



rh q d£ 



d 1 (^;A,B,C)' 



R\A,B,C,s) = 2 / ^ 



(Z;A,B,C)' 



T(x,y,z;A,B,C,s)=2 



n(2-p),b n £ 



1/;: 



H l I (x,y,z;A,B,C,s)=2 



d^-A,B,C) 



logm g g!£ 



(31) 

(32) 
— ,(33) 



di(f;A,B,C)' 



H^(x,y,z;D,A,B,C,q,s)=2 J 



S(x,y,z;D,A,B,C,s) = 2 



m q t 2 d£ 



d 2 {i;D,A,B,C)' 



(34) 



(35) 



d 2 (£;D,A,B,C) 



(36) 



where 



- 2 ,2 
m — E, 



+ 



+ 



l + (a 2 - 1)£ 2 1 + (6 2 - 1)£ 2 1 + (c 2 - 1)5 



(37) 



and 



di(C;A,B,C) = 

V(l + (A- 1)C 2 )(1 + (5 - 1)? 2 )(1 + (C - 1)£ 2 ), (38) 

d 2 (£; D, A, B, C) = (1 + (D- 1)£ 2 ) A, B, C). (39) 

Total mass of the galaxy represented by our new triaxial 
model is easily found to be 



Mtg = 2a 



3-7 



pRlnb n{v ~ 3) V 



l/r. 



.(40) 



To "normalize" our model to have the total mass equal to 
unity, the expressions for potential and forces need to be 
divided by Mtg- 

In the spherical limit, a = b = c = 1, the expressions of 
Terzic & Graham (2005) are recovered (see Appendix A). 



2.3 Special Case r% = 0: Triaxial Generalization of 
Prugniel— Simien Profile 

In the limit of the break radius rb going to zero, our two- 
regime model reduces to a single regime Sersic profile. The 
corresponding equations for mass density, potential and 
forces for this triaxial profile are readily obtained after set- 
ting rb = in the equations (f2}, (|ll|l and (|26[) . respectively, 
and are given by 



p{m) = p b p y-j^- J e 



(41) 



$(x,y,z) = -anpRlbZ (p 2> T(x, y, z; a 2 , b 2 , c 2 , 00), (42) 
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F qi {x,y,z) = -apqiS(x,y,z;ki,a ,b ,c ,00) 



(43) 



Total mass for the Prugniel-Simien model is readily ob- 
tained by substituting r b = into (l40t : 



M PS = 2apRlnb n(p - 3) r [n(3 - p)] 



(44) 



2.4 Special Case rt — 00: Triaxial Power-Law Core 

In the limit of the break radius rt tending to infinity, our 
two-regime model reduces to a single power-law core profile. 
The equations for mass density, potential and forces for the 
3D power-law core are given by 



p(m) = Pb\-^- 
\ m 



(45) 



'ffOE,y, 2 ;a 2 ,6 2 ,c 2 ,2- 7 ), (46) 



-arlq l H 2 {x, y, z; k 2 ,a 2 , b 2 , c, -7), (47) 



above expressions are given in terms of elementary func- 
tions and Carlson integrals ( |de Zeeuw fe P fcrmigc r~l988| 
IPoon fe Merritt 200 ljl . Single power-law model is scale-free, 
which means that orbital properties simply scale with ra- 
dius. This property leads to the separation of radial and 
angular dependence in spherical coordinates, thus allowing 
for efficient series approximation, such as double Fourier or 
Fourier- Legendre {Tcrzic 2002} ■ 

Total mass for the power-law core is readily found to be 



2 « 3 
M = r 

Ay± o max ! 

3-7 

where r max is the outer radius of the core. 



(48) 



2.5 Non-Dimensionalizing the Physical Problem 

In order to make the numerical problem dimensionless, we 
adopt a convention G = Mtot = Mtg — Mps — Md = 1. 
This means that, since 



G = 6.672 x 10 _11 [kg] _1 [m] 3 [sec] -2 
= [lO 11 M ]" 1 [kpc] 3 [1.491 x 10 6 yif 

the unit of time is 



1.491 x 10 V /MQ 1011 ' 



M \ lkpc J ' 



(49) 



(50) 



where j3 is the length-scale of the model. The total stellar 
mass can be easily computed from the best-fit parameters 
of three models, via equations (|40)l . (|44)l and (|C2[1 . If the 
observed effective half-light radius is taken to be the length- 
scale of the model j3 (see Tabled, respective time-scales are 
easily computed using expression (|50[) . Each unit of length 
corresponds to 1 kpc. 



3 NUMERICAL IMPLEMENTATION 

In this Section we outline the major features of the nu- 
merical suite which we designed to optimize orbit integra- 
tions in the new triaxial potential - triaxial generalization 



of a 5-parameter Terzic-Graham model, as well as its 3- 
parameter special case (when central core is not depleted) 
triaxial Prugniel-Simien model. 

Our goal is to design a state-of-the-art, user-friendly 
numerical package which would allow for easy setup and 
efficient execution of numerical integrations of orbits in 
the new model. Our suite of numerical routines is thor- 
oughly tested and fully optimized. It is also parallelized us- 
ing MPI, with its efficiency scaling roughly linearly with the 
number of processors used. We make this numerical suite 
freely available for downloading from our group's webpage: 
http : //www.nicadd.niu. edu/research/beams/TS_2007/ 



3.1 Outline of the Numerical Suite 

The numerical suite utilizes a number of tested and 
optimized non-proprietary numerical packages writ- 
ten in FORTRAN 77. The numerical algorithm used 
for integrating orbits is the explicit imbedded (7,8) 
pair of Dormand and Prince to compute orbital points 
( |Hairer, Norsett fe Wanner 1993[ ). Complete and incom- 
plete gamma functions are computed by routines from 
SLATEC, created by American National Laboratories 
(available online at http://www.netlib.org/slatec). 
Evaluation of quadratures, both with finite and 
infinite ranges, is done by QUADPACK routines 
(Piessens, de Doncker-Kapenga & Uberhuber 19831, (also 



available online at http://www.netlib.org/quadpack). 

The driver program is written in FORTRAN 90. 

Flow-chart outlining the numerical suite is shown in 
Fig.H 



3.2 Modular Design 

We designed the integrator core of the numerical package to 
be as simple as possible while allowing easy addition of mul- 
tiple models and analysis routines in a parallelized integrator 
system. The data flow of the system is outlined in Fig. H A 
parameter settings file is defined by the user, which defines 
settings for the integrator, the initial conditions generator, 
the potential and forces models, and the analysis routines. 
This file is first processed by the initial conditions generator 
to produce an initial conditions file. The integrator core then 
uses both the parameter settings file and the initial condi- 
tions file to begin integrating the orbits. For each orbit, the 
initial conditions are read directly from the initial condi- 
tions file, and the integrator core obtains the forces from 
the models using the standard "potential and forces" inter- 
face in order to perform the total orbital integration. When 
the integration is complete, orbital data is passed on to the 
analysis routines, which are responsible for data processing 
and producing the output file for the simulation. 

Every effort is made to separate the user-editable code 
from the system-level code. This allows the system to be 
modified for use with any model and analysis routines with- 
out requiring changes to the source files containing the in- 
tegrator or parallelization code. The user can supply any 
model to the integrator using the standard "potential and 
forces" interface. This is done by simply filling in the sup- 
plied functions get_acceleration and get_potential with 
user-specific code. Supplying analysis routines such as an 
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Figure 2. Flow-chart outlining the numerical suite. 



FFT routine or an energy-conservation measurement is sim- 
ilarly trivial, by simply adding code to perform the operation 
in the analyze_orbit routine. More advanced analyses such 
as tracking Lyapunov exponents or tangent dynamics are ob- 
tained through a similar interface which allows an arbitrary 
number of extra integration variables to be tracked along 
with the main orbit without changes to the main integrator 
source code. 

We created a parallel implementation of this integra- 
tor using MPI by using a block-decomposition technique 
to assign nearly equal amounts of orbits to each processor. 
Each processor computes its orbits one at a time, and writes 
its output to its own local output file. Once all processors 
have computed their sets of orbits, the main output file is 
constructed by concatenating the individual output files to- 
gether. The final output file contains the orbits in exactly 
the same order as they were arranged in the initial condi- 
tions file because a block-decomposition assignment of orbits 
is used. The parallel implementation is completely transpar- 
ent to the user-editable code in the analysis and model def- 
inition routines, because the parallel decomposition of the 
problem is done in groups of orbits in the same manner as 
the analysis decomposition. 



3.3 Parallel Implementation 

Our parallel decomposition of the orbital integrator requires 
very little communication between the computation nodes, 
which is why we estimate that the speedup derived from 
parallel computation scales roughly linearly with the num- 
ber of processors used. However, this nearly optimal speedup 
is limited by two main factors. First, the number of orbits to 
integrate must exceed the number of processors. This is re- 
quired in order to allow each processor to contribute toward 
the overall solution. Second, the length of time to integrate 



each orbit must be approximately equal, because the algo- 
rithm is not finished until the slowest processor has com- 
pleted its equal sized block of orbit integrations. 

In common use of this integrator, however, these factors 
are not expected to play a large role in the actual perfor- 
mance of the simulation. The first problem is averted by the 
fact that integration batches are commonly grouped in larger 
numbers than are typical for a multiple-processor cluster. 
Common clusters have processor counts in the hundreds, 
while it is often desirable to run integration groups num- 
bering from hundreds to thousands of stars in each energy 
level. The second problem is avoided because it is common 
practice to group initial conditions in batches according to 
energy level. Since the dynamical time is roughly equal for 
all stars of the same energy level, we can say that all in- 
tegrations of the same simulation time length will take ap- 
proximately the same amount of real time to simulate. One 
exception to this is the fact that centrophilic orbits, ones 
with no angular momentum, typically take longer to inte- 
grate: the integrator is forced to take increasingly smaller 
steps near the very center in order to attain prescribed ac- 
curacy when orbits' velocities are high. This is why in our 
numerical suite generates initial conditions for centrophilic 
orbits (vanishing angular momentum) separately from the 
centrophobic ones (see Appendix B). Therefore, if these two 
types of orbits are integrated separately, the roughly linear 
scaling of simulation time with the number of processors is 
preserved. 



4 MODEL COMPARISON 

The Dehnen model and simple power-law core models 
have been quite popular as qualitative, "toy" models which 
capture general features of stellar dynamics near galactic 
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centers (|Merritt fc Fridman 19961 IValluri fc Merritt 19981 
IMerritt fc Valluri 1999"1lPoon fc Merritt 20011 ITerzic 2002p . 
By devising parametric models which provide excellent fits 
to the full range of the observed data, not only the central 
region, we are now in position to develop quantitative models 
for individual galaxies. 

Here we propose the 3-parameter triaxial generaliza- 
tion of Prugniel-Simien model as an alternative to the 3- 
parameter Dehnen model for modeling dynamics of ellipti- 
cal galaxies without depleted cores. For the galaxies with 
depleted cores, we propose a triaxial generalization of the 
5-parameter Terzic-Graham. We first reiterate the findings 
of Terzic & Graham (2005) that these alternatives provide 
much better fits to deprojected observed light-profiles of re- 
alistic galaxies using two representative examples from their 
set of eight elliptical galaxies and the bulge of the Milky 
Way. We then demonstrate that using these improved tri- 
axial models in dynamical simulations will not incur ap- 
preciable additional computational cost over the traditional 
Dehnen model. It is our belief that having a consider- 
ably more accurate model of the observed physical system 
more than justifies the modest increase in computational 
expenses. 

4.1 Fitting Realistic Deprojected Light-Profiles 

Terzic & Graham (2005) fitted several of the most pop- 
ular density models: 2-parameter Jaffe (1980) and Hern- 
quist (1990); 3-parameter Dehnen and Prugniel-Simien; and 
their own new 5-parameter hybrid model to deprojected 
light profiles of several elliptical galaxies. For their sample 
of luminosity-density profiles taken from real galaxies the 
2-parameter Jaffe and Hernquist double power-law mod- 
els systematically fail to fit the entire extent of observed 
galaxy profiles. The 3-parameter Dehnen model does bet- 
ter, providing a good match to density profiles of large 
early- type galaxies with Sersic index around 4 or greater, 
but it is inadequate for galaxies with low Sersic indices, 
which include dwarf ellipticals and most bulges in disk galax- 
ies. The 3-parameter Prugniel-Simien model, however, ac- 
curately matches the density profiles of both luminous ellip- 
tical galaxies with n > 4 and faint elliptical galaxies with 
n < 4, including exponential (n — 1) profiles. For galaxies 
with partially depleted cores, Terzic & Graham (2005) de- 
veloped a 5-parameter, two-regime model, capable to match 
simultaneously both the nuclear and global stellar distribu- 
tion. Terzic & Graham (2005) arrive to these conclusions 
after fitting a set of eight elliptical galaxies of various sizes 
and light-profile shapes. 

In the present study, we focus on the bulge in the cen- 
ter of the Milky Way and the two galaxies from the sample 
used in Terzic & Graham (2005): an elliptical galaxy NGC 
1379, and an elliptical galaxy with partially depleted core 
- NGC 3348. The Milky Way's bulge and the first galaxy, 
NGC 1379, with Sersic index of n = 2.0 (see Table 2 of 
Terzic & Graham 2005), were chosen because it illustrates 
that Dehnen is inadequate to describe galaxies with low 
Sersic indices, and that Prugniel-Simien model is a good 
alternative. The second galaxy, NGC 3348, is an example 
of a galaxy with partially depleted cores for which only 5- 
parameter Terzic-Graham model provides good fits in den- 
sity space. 



The profiles have been calibrated in units of Lq pc 
using the distance moduli provided in Tonry et al. (2001) 
for NGC 1379 and the distance of 41 Mpc for NGC 3348. 
We used solar absolute magnitudes of Mk = 3.33, Mb = 
5.47 and M R = 4.28 (Cox 2000). A Hubble constant H = 
75km s _1 Mpc -1 is used. 

In the top two rows of Fig. [3] the 3-parameter Dehnen 
and Prugniel-Simien models were fitted to the deprojected 
major-axis K-band light-profile of the Milky Way's bulge 
(first row) and B-band light-profile of NGC 1379 (second 
row). It is clear that Dehnen model is morphologically in- 
adequate to provide a good fit for these types of galaxies 
which have a low value of Sersic index n. Prugniel-Simien 
model does better, especially fitting the outer curvature 
of the profile, as evidenced by the smaller relative error 
5{r). In the bottom row of Fig. [3] the 3-parameter Dehnen 
and Prugniel-Simien models, along with the 5-parameter 
Terzic-Graham were fitted to the deprojected major-axis 
i?-band light-profile of a core elliptical galaxy NGC 3348. 
Both Dehnen and Prugniel-Simien models fail to account for 
the inner partially depleted core. The 5-parameter Terzic- 
Graham model does quite well and is the only density model 
that can fit galaxies possessing depleted cores. 

Fig. [4] shows the total enclosed stellar mass, scaled po- 
tential and forces for the models from Fig. [3j using the best- 
fitting parameters reported in Table [2] 

4.2 Computational Efficiency 

The quadratures involved in computing forces for the 
Dehnen model are similar to those used for Prugniel-Simien 
and Terzic-Graham models (compare the equation ()C6[) to 
equations (|25l) - (|26l ) and (1431) '). It is therefore reasonable to 
expect that alternative implementations of the numerical 
evaluation of these quadratures will be equally reflected on 
the computational speeds. In other words, the relative com- 
putational efficiency outlined later in the Section should be 
general. 

It is possible to considerably improve computational ef- 
ficiency of force computations by exploiting simplifications 
arising for special cases of parameters (e.g. integer values of 
the central cusp 70 for the Dehnen model) or galaxy shapes 
(e.g. spherical case a = b = c = 1). These simplifications 
were exploited in the earlier, qualitative studies featuring 
the Dehnen model, in which 70 = represented the model 
with no central cusp, 70 = 1 "weak" and 7b = 2 "strong" 
central cusp (Friedman & Merritt 1996, Siopis 1999). How- 
ever, for qualitative studies such as the one we are proposing 
here, it is not expected for deprojected galaxy profiles to be 
best fitted by parameter values which take on these fortu- 
itous values, which is why these simplifications arising from 
special cases are not explored here. 

4.3 Comparing Speed of Orbital Integration 

We now compare numerical efficiency of orbital integrations 
of the two new models - Terzic-Graham and its special case 
Prugniel-Simien - with the traditional Dehnen model for a 
typical galaxy. (Expressions for density, potential and forces 
for the Dehnen model are given in Appendix C.) We compare 
integration times for equivalent sets of orbits in the three 
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Figure 3. Top row: parametric fits (along with the relative error shown underneath) of Dehnen model (left), Prugniel-Simien model 
(right) to the deprojected major-axis, fC-band light profiles of the central bulge in the Milky Way (fitted out to R = 9 degrees; data 
from Kent, Dame & Fazio 1991). Middle row: parametric fits (along with the relative error shown underneath) of Dehnen model (left), 
Prugniel-Simien model (right) to the deprojected major-axis, B-band light profiles of NGC 1379 elliptical galaxy (data from Caon, 
Capaccioli & D'Onofrio 1993). Bottom row: parametric fits (and relative errors) of Dehnen model (left), Prugniel-Simien model (middle) 
and Terzic-Graham (right) to deprojected major-axis, ij-band light profiles of NGC 3348 elliptical galaxy (data from Trujillo et al. 2004). 
The luminosity-density profile u(r) = p(r) / (M / L) is in units of Lq pc -3 . The model parameters are given in Table [2] 
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Table 1. For the bulge at the center of the Milky Way, NGC 1379 and the core galaxy NGC 3348: total stellar mass, unit of model 
(simulation) time in physical time and the Hubble time (1.37 X 10 10 years) in terms of model (simulation) time, as computed from the 
parameters for the three models via equation J50D - Wc adopt the observed half-light radius computed from Scrsic and core-Sersic fits 
to projected light-profiles, reported in Table 2 of Terzic & Graham (2005), as the length-scale /3. Mass-to-light ratios are taken to be 
M/L = 0.66 for the iv-band, M/L = 3 for the R-band, and M/L = 5.3 for the B-band (Worthey 1994). 



Gal. Id. 


P 


Model 


Total stellar mass 


Unit of time 


^Hubble 




[arcsec] 




[%] 


[years] 


[model units] 


Milky Way 


1.65 x 10 4 


Dehnen 


1.00 x 10 u 


4.60 x 10 b 


2.98 X 10 3 






Prugniel-Simien 


5.34 x 10 8 


6.30 x 10 6 


2.17 X 10 3 


NGC 1379 


24.3 


Dehnen 


5.93 x 10 iu 


4.50 x 10 b 


3.04 x 10 3 






Prugniel-Simien 


5.02 x 10 11 


4.90 x 10 6 


2.78 x 10 3 


NGC 3348 


21.4 


Dehnen 


1.86 x 10 11 


6.11 x 10 b 


2.24 x 10 3 






Prugniel-Simien 


1.76 x 10 11 


6.29 x 10 6 


2.18 x 10 3 






Terzic-Graham 


1.72 x 10 11 


6.37 x 10 6 


2.15 x 10 3 



Table 2. Best-fitting parameters from the three density models over the entire range of the deprojected profile. The units are Lq pc 3 
for p and arcseconds for R e , r a and r b . 



Gal. Id. Band — Dehnen — — Prugniel-Simien — — Terzic-Graham 







r a 


log p(r a ) 


ID 


Rc 


logp(Ac) 


n 


Re 


n 7 n 


log p(r b ) 


Milky Way bulge 


K 


1.20 x 10 4 


-7.31 


0.00 


1.57 x 10 4 


-7.52 


1.08 








NGC 1379 


B 


11.1 


-0.49 


0.00 


24.7 


-1.31 


2.10 








NGC 3348 


R 


6.40 


0.14 


0.71 


13.2 


-0.63 


2.10 


20.2 


3.6 0.44 0.37 


2.15 



models obtained after parametric fitting in density space to 
the deprojected light-profile of NGC 3348 (Table |3). 

In comparing the computational speeds of various mod- 
els against one another, it is necessary to determine a setting 
which provides a fair comparison. On the surface, it seems 
that comparing run times for orbital integrations of equal 
integration time would be a reliable setting. In fact, this is 
the case when using a constant time-step ordinary differen- 
tial equation (ODE) solver, because it guarantees that the 
get_acceleration routine will be called an identical number 
of times for each model. 

However, restricting the orbital integrator to con- 
stant time-step solvers is not desirable because of their 
poor performance when compared to adaptive-step solvers. 
Adaptive-step solvers adjust the size of the integration step 
so as to satisfy the predefined error tolerance. Simulation 
time for integrations using adaptive-step ODE solver only 
obtains physical significance when scaled to dynamical time 
of the orbit at a certain energy level. In this setting, we 
define the dynamical time to be the average period of an 
oscillation of a star around the origin of the galaxy. 

For the NGC 3348 (parameters are given in Table[2]), we 
integrated 600 box orbits and 600 tube orbits (see Appendix 
B) at 50 logarithmically spaced energy levels in the range 
E{lQ- 2 f3) < E{r) < £7(5/3) (where E(r) = $(r,0,0) denotes 
a isoenergy surface pierced by the x-axis at x = r), for 
roughly 3 Hubble times (length-scale /3 and Hubble times are 
given in Table [!}. We used two axis ratios, that of spherical 
model a = b = c = 1, and that of a maximally triaxial 
((a 2 - 6 2 )/(a 2 - c 2 ) = 0.5) a : b : c = 1 : 0.79 : 0.5. 

Figure[5]shows the total simulation time. The Prugniel- 
Simien is consistently about 40 — 80% slower than the 
Dehnen model, while Terzic-Graham, on average, is at least 
as fast as the Dehnen model. There is a "bump" in the 
line corresponding to the Terzic-Graham model around the 
break radius. This is due to the fact that the adaptive-step 
ODE integrator samples the region around the transition 



more densely, which results in increase in the total simula- 
tion times. 

Increasing the steepness of the negative central den- 
sity slope (7_d for Dehnen, 7 for Terzic-Graham and p 
for Prugniel-Simien) of the model renders orbital integra- 
tion more computationally intensive, because of the rapid 
changes in the forces. For NGC 3348, the central density 
slopes of the Dehnen (70 = 0.71) is nearly identical to that 
of the Prugniel-Simien (p — 0.72), but somewhat larger than 
that of Terzic-Graham (7 = 0.44). This means that the com- 
parison of execution speeds for orbits in potentials modeling 
NGC 3348 will accurately reflect efficiency of force evalua- 
tions in the Prugniel-Simien model relative to that of the 
Dehnen model, while the Terzic-Graham will have a slight 
"advantage", having the shallower central slope. This is in- 
deed what we see in Fig. \5\ the 5-parameter Terzic-Graham 
model outperforms both 3-parameter models - the tradi- 
tional Dehnen and the Prugniel-Simien - mostly because of 
its shallower central cusp. 

The full 5-parameter version of the Terzic-Graham 
model should only be used to model galaxies with depleted 
cores, which means that it would always have shallower in- 
ner slopes than their Dehnen or Prugniel-Simien counter- 
parts, and therefore be computationally very competitive 
with them, often even outperforming them (the exact per- 
formance depending on the level of depletion of the core). 

The special case of the 5-parameter Terzic-Graham, the 
3-parameter Prugniel-Simien, should only be used to model 
galaxies without depleted cores. In those cases, the com- 
putational price to pay over the traditional Dehnen model 
is modest - about 40 — 80% over the meaninful range of 
energies. We believe that having a quantitative model ap- 
preciably more faithful to the physical problem is well worth 
the price. 
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Figure 4. Top column: mass enclosed within radius r in terms of Mq (top), scaled potential (middle) and force (bottom), along the 
x-axis, for a spherical Dehnen model (dashed lines) and Prugniel-Simien model (dotted lines) for the Milky Way bulge. Middle column: 
mass enclosed within radius r in terms of Mq (top), scaled potential (middle) and force (bottom), along the x-axis, for a spherical Dehnen 
model (dashed lines) and Prugnicl-Simicn model (dotted lines) for NGC 1379 elliptical galaxy. Right column: mass enclosed within radius 
r in terms of Mq (top), scaled potential (middle) and force (bottom), along the x-axis, for a spherical Dehnen model (dashed lines), 
Prugnicl-Simicn model (dotted lines) and Terzic-Graham (solid lines) for NGC 3348 elliptical galaxy. The model parameters are given 
in Table [2] The axis ratios are 1:1:1. 



5 SUMMARY 

Motivated by the findings of earlier studies (Graham et 
al. 2002, 2003; Terzic & Graham 2005) which demonstrate 
that the traditional power-law density models, such as 3- 
parameter Dehnen, are inadequate to describe density pro- 
files whose slopes continuously vary as a function of ra- 
dius - as is the case for most elliptical galaxies and bulges 
- we introduce a new family of realistic triaxial density- 
potential-force profiles for stellar systems and dark matter 
halos. The new family is a 5-parameter, two-regime model 
of inner power-law density core and the outer deprojected 
Sersic profile, with the boundary between the two regimes 
being the break radius 7V For extreme values of the break 
radius the model reduces to the 3-parameter deprojected 
Sersic R 1 '" profile (rt — ► oo) and the 2-parameter power- 
law core (rt — > 0). The full 5-parameter model is invoked 
only when modeling elliptical galaxies with depleted cores; 
in other cases, which comprise the majority of observed el- 
liptical galaxies and bulges, the 3-parameter triaxial depro- 
jected Sersic model suffices. 



In Section [21 we derived the potential and forces for the 
fully triaxial density model first presented in its spherical 
form by Terzic & Graham (2005). Potential and forces are 
expressed in terms of easy-to-evaluate quadratures. In Sec- 
tion|3] we developed an optimized suite of numerical routines 
designed for use in orbital integrations. 

The goal of this study was to develop and offer to the 
community an alternative to the traditional Dehnen model 
which is substantially more faithful to the physical prob- 
lem (as measured by the quality of fits to deprojected light 
profiles) at an almost marginal computational expense. 
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Figure 5. Total simulation time for 600 orbits at varying energy 
levels (denoted by the maximum radial excursion) in potentials 
approximating NGC 3348 galaxy, integrated for the equivalent of 
3 Hubble times. Top row: box orbits (left) and tube orbits (right) 
in the spherical geometry (axis ratio 1:1:1). Bottom row: 
box orbits (left) and tube orbits (right) in the maximally triaxial 
geometry (axis ratio 1 : 0.79 : 0.5). Crosses represent Dehnen 
model, triangles Prugniel-Simien and circles Terzic-Graham. 
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APPENDIX A: SPHERICAL LIMIT 

In the spherical limit a — b = c — 1, m = r and m(r) = 
f(r) = r/s/r + 1 and all quadratures above have analytic 
solutions given in terms of elementary and special functions: 



H'(x, y, z; 1, 1, 1, q, s) = r" J (r + 1)" 



dT 



l + 2q 

H(x,y,z;l,l,l,q) = tffay, z; 1, 1, l,q, 0) 



(s + !)-?-«, (Al) 
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The result for T(x,y, z; 1, 1, 1, s) is derived after integration 
by parts, tj, is easily found from = r/y/r + 1. 

After some straightforward yet tedious algebra, the 
equations for mass density, potential and forces reduce to 
those of Terzic & Graham (2005), their equations (6), (8) 
and (14), respectively. 



APPENDIX B: 
GENERATOR 



INITIAL CONDITIONS 



Sampling of the initial condition space is done by em- 
ploying a two-fold two-dimensional start space, following 
Schwarzschild (1993): the stationary start space which con- 
tains initial conditions starting from equipotential surfaces 
with zero velocities; and the principal-plane start space 
which consists of radially stratified initial conditions which 
pierce one of the three principal planes with the velocity 
vector normal to the plane. These start spaces are designed 
to sample different types of orbits arising in triaxial poten- 
tials: stationary start space picks up orbits which have zero- 
velocity turning points, such as boxes and resonant boxlets, 
while the principal-plane start space selects mostly tube or- 
bits (|Schwarzschild 19931 ITerzic 2002p . The two start spaces 
are shown in Fig. IB1I 

The spherical angles of the stationary start space loca- 
tions are chosen as in Schwarzschild (1993), with the addi- 
tion that the number of points can be chosen using a pa- 
rameter statJJ, so that the total number of points on this 
equipotential surface will be 3 statJJ 2 . The zero- velocity ra- 
dius of this potential in the given direction is then found, 
resulting in an exact coordinate for the given initial condi- 
tion. 

We calculated zero- velocity radii using a secant-method 
search to find the root of $(r, 6,<j>) — E = 0, given the di- 
rection angles 6 and (f>. This method results in very fast 
convergence to the solution, even on the unbounded interval 
< r < 00. Proper convergence of this method does as- 
sume, however, that the <I>(r, 0, (f>) function is monotonically 
increasing with radius r, which is a condition that is met by 
the three potential functions considered in this study. 

Initial conditions in the principal-plane start space are 
calculated in a similar manner to Schwarzschild (1993) as 
well. The parameters plane_N_theta and plane_N_r deter- 
mine the number of points sampled. First, the angle 6 is 
subdivided into plane_N_theta regions from to 7r/2. The 
centers of each of these regions are chosen as direction vec- 
tors within the plane for choosing initial conditions. Along 
each direction vector, the zero- velocity radius rmax is com- 
puted as described above. Next, an inner cutoff radius is 
determined by using a fraction plane_f rac_rmin such that 
rmin = plane_f racjrmin X rmax. This inner cutoff fraction 
is a parameter that can be chosen from to 1 in order to re- 
duce the number of orbits that are duplicated at small radii. 
The initial conditions are then chosen by subdividing the re- 
gion rmin to rmax into planeJLr sections, and selecting the 
center of each section as the location of the initial condi- 
tion. The velocity is then chosen normal to the plane such 
that the total energy is equal to the desired energy level. 
This is done for each of the 3 principal planes, resulting in 
3 x plane_N_theta x plane_N_r initial conditions generated 
in this start space. 



Triaxial density-potential-force profiles 




Figure Bl. Stationary and principal-plane start spaces. 



APPENDIX C: DEHNEN MODEL 



1 m-"< D {r a + my D - A i 2 



Dehnen model (|Dehnen 19931 ITremaine et al. 1994]) is a J d 2 (£,; k 2 ; a 2 ,b 2 , c 2 ) 

double-power-law, 3-parameter density model, defined by its 

steepness of the inner cusp 70, break radius r a which marks where fci = a, k 2 = b, k 3 = c and gi = x, q 2 = y, qs = 
the transition between the two power-law regimes, and the 
density at the break radius p(r a ): 



p(m) = 2*-">p(r„) (— J (r a + m)^-\ (CI) 
The total mass is easily found to be 

M D = ^Lrl2^°p{r a ). (C2) 

3 — 7d 

The corresponding potential is given by 

00 

\ GMd f P(rh) 
${x,y,z) = -— — / — v 2 ' 2 dr, C3) 

2(2-7rj)r a J di(r; a 2 , b 2 , c 2 ) 


where 

, / ffl \ 2 -TD 

P(m) = 1 - (3 - 70) 



r a +rn 
1J1 



(2- 7 d)(— ^-)^ 7D (C4) 

This equation is, again, transformed into a proper integral 
via transformation £ = (r + 1)~ 1,/2 , to obtain 



a/ \ GMd /" P(m) 

$(x,y,z) = — / v y d£. (C5) 



Similarly, the forces are found to be, after the change of 
variables, 

F qi (x,y,z) = -GMd(3 - jD)r a qi 



